{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# An Example of using sklearn Pipeline with matminer\n",
    "\n",
    "This goes over the steps to build a model using sklearn Pipeline and matminer. Look at the intro_predicting_bulk_modulus notebook for more details about matminer and the featurizers used here.\n",
    "\n",
    "This notebook was last updated 11/15/18 for version 0.4.5 of matminer.\n",
    "\n",
    "**Note that in order to get the in-line plotting to work, you might need to start Jupyter notebook with a higher data rate, e.g., ``jupyter notebook --NotebookApp.iopub_data_rate_limit=1.0e10``. We recommend you do this before starting.**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Why use Pipeline?\n",
    "\n",
    "Pre-processing and featurizing materials data can be viewed as a series of transformations on the data, going from the initially loaded state to training ready. Pipelines are a tool for encapsulating this process in a way that enables easy replication/repeatability, presents a simple model of data transformation, and helps to avoid errant changes to the data. Pipelines chain together transformations into a single transformation. They can also be used to build end end-to-end methods for preprocessing/training/validating a model, by optionally putting an estimator at the end of the pipeline. See the [scikit-learn Pipeline documentation](http://scikit-learn.org/stable/modules/generated/sklearn.pipeline.Pipeline.html) for details."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load sklearn modules\n",
    "from sklearn.pipeline import FeatureUnion, Pipeline\n",
    "from sklearn.base import TransformerMixin, BaseEstimator\n",
    "\n",
    "from sklearn.linear_model import LinearRegression\n",
    "from sklearn.ensemble import RandomForestRegressor\n",
    "from sklearn.svm import SVR, LinearSVR\n",
    "\n",
    "from sklearn.decomposition import PCA, NMF\n",
    "from sklearn.feature_selection import SelectKBest, chi2\n",
    "from sklearn.preprocessing import StandardScaler, MinMaxScaler\n",
    "\n",
    "from sklearn.metrics import mean_squared_error\n",
    "from sklearn.model_selection import RepeatedKFold, cross_val_score, cross_val_predict, train_test_split, GridSearchCV, RandomizedSearchCV\n",
    "\n",
    "import numpy as np\n",
    "from pandas import DataFrame\n",
    "from scipy.stats import randint as sp_randint\n",
    "\n",
    "# Load featurizers and conversion functions\n",
    "from matminer.featurizers.composition import ElementProperty, OxidationStates\n",
    "from matminer.featurizers.structure import DensityFeatures\n",
    "from matminer.featurizers.conversions import CompositionToOxidComposition, StrToComposition"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Loading the Dataset\n",
    "Matminer comes pre-loaded with several example data sets you can use. Below, we'll load a data set of computed elastic properties of materials which is sourced from the paper:  \"Charting the complete elastic properties of inorganic crystalline compounds\", M. de Jong *et al.*, Sci. Data. 2 (2015) 150009."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matminer.datasets.convenience_loaders import load_elastic_tensor\n",
    "df = load_elastic_tensor() # loads dataset in a pandas DataFrame \n",
    "unwanted_columns = [\"volume\", \"nsites\", \"compliance_tensor\", \"elastic_tensor\", \n",
    "                    \"elastic_tensor_original\", \"K_Voigt\", \"G_Voigt\", \"K_Reuss\", \"G_Reuss\"]\n",
    "df = df.drop(unwanted_columns, axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# seperate out values to be estimated\n",
    "y = df['K_VRH'].values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data Preprocessing\n",
    "The conversion functions in matminer need to be run before the pipeline as a data preprocessing step."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "8100ec8295fd4ef49f3ec36bc796d25c",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, description='StrToComposition', max=1181), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "f72c6d834103430e81fb0d48376bc81e",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, description='CompositionToOxidComposition', max=1181), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "df = StrToComposition().featurize_dataframe(df, \"formula\")\n",
    "df = CompositionToOxidComposition().featurize_dataframe(df, \"composition\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Helper Functions\n",
    "The matminer library uses pandas DataFrames, where sklearn.pipeline mainly looks at things as numpy arrays, so helper methods are needed to seperate out columns from the DataFrame for pipeline. To be used in pipeline they need to be transformers, meaning they implement a transform method. (A fit method that does nothing is also needed)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matminer.utils.pipeline import DropExcluded, ItemSelector"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Making Feature Union Pipeline for Featurizers\n",
    "This creates a pipeline that transforms preprocessed data to featurized data usable in sklearn. It can be used to transform data on its own or as part of another pipeline. It is possible to cache values in the pipeline so that this is only done once.\n",
    "\n",
    "This Feature Union pipeline has three parts, ``drop`` which drops unwanted columns, ``density`` which adds density features, and ``oxidation`` which adds oxidation state features. These are combined by ``FeatureUnion`` to create the final dataset. The ``drop`` transform acts as an identity+filter, passing through the original data minus the unwanted columns. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# columns to remove before regression\n",
    "excluded = [\"G_VRH\", \"K_VRH\", \"elastic_anisotropy\", \"formula\", \"material_id\", \n",
    "            \"poisson_ratio\", \"structure\", \"composition\", \"composition_oxid\"]\n",
    "\n",
    "# featurization transformations\n",
    "featurizer = FeatureUnion(\n",
    "    transformer_list=[\n",
    "        ('drop', DropExcluded(excluded)),\n",
    "        ('density', Pipeline([\n",
    "            ('select', ItemSelector(\"structure\")),\n",
    "            ('density_feat', DensityFeatures())\n",
    "        ])),\n",
    "        ('element', Pipeline([\n",
    "            ('select', ItemSelector(\"composition\")),\n",
    "            ('oxidation_feat', ElementProperty.from_preset(preset_name=\"magpie\"))\n",
    "        ])),\n",
    "        ('oxidation', Pipeline([\n",
    "            ('select', ItemSelector(\"composition_oxid\")),\n",
    "            ('oxidation_feat', OxidationStates())\n",
    "        ])),\n",
    "    ]\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Making a Regression Pipeline\n",
    "This is a simple pipeline combining the featurizer transformer pipeline with a linear regression estimator."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "e90201712c484a6f97830870dca64003",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, description='DensityFeatures', max=1181), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "0a4654cd84fa4ee488e88efa071980a6",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, description='ElementProperty', max=1181), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "bcb0b0bcde8544738b758b469c1a1b57",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, description='OxidationStates', max=1181), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/daniel/.conda/envs/matminer/lib/python3.6/site-packages/sklearn/linear_model/base.py:485: RuntimeWarning:\n",
      "\n",
      "internal gelsd driver lwork query error, required iwork dimension not returned. This is likely the result of LAPACK bug 0038, fixed in LAPACK 3.2.2 (released July 21, 2010). Falling back to 'gelss' driver.\n",
      "\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "b4235a9f0aae45baa3729e98f6cc3e55",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, description='DensityFeatures', max=1181), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "0e6cd266e8e84218995780c9bb2b3b07",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, description='ElementProperty', max=1181), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "84ca9dfc400f4ec5861fed22d2517701",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, description='OxidationStates', max=1181), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "training R2 = 0.928\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "77a7fafb457f45c6b42c4a333608bd10",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, description='DensityFeatures', max=1181), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "f1f0112578224b318db718e1a205d638",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, description='ElementProperty', max=1181), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "331186822f8946eb935fa0a161aad7ea",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, description='OxidationStates', max=1181), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "training RMSE = 19.595\n"
     ]
    }
   ],
   "source": [
    "# make the pipeline\n",
    "pipeline = Pipeline([\n",
    "    ('featurize', featurizer),\n",
    "    ('regress', LinearRegression()),\n",
    "])\n",
    "\n",
    "pipeline.fit(df, y)\n",
    "\n",
    "# get fit statistics\n",
    "print('training R2 = ' + str(round(pipeline.score(df, y), 3)))\n",
    "print('training RMSE = %.3f' % np.sqrt(mean_squared_error(y_true=y, y_pred=pipeline.predict(df))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Making a Random Forest Pipeline\n",
    "This is the same, but with a random forest regression instead. The only line changed is the one defining ``regress`` in the pipeline."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "3249dac722f34ccab79ceb3b1b000682",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, description='DensityFeatures', max=1181), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "8e74864eeccc41a09511c74c72d2b085",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, description='ElementProperty', max=1181), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "dbb8ac9bbf3740f6831e779e0d6bc697",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, description='OxidationStates', max=1181), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "e85a9198729644888fe73a10d677732e",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, description='DensityFeatures', max=1181), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "9f0953d21c8349e493aeb4b3da13e7ee",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, description='ElementProperty', max=1181), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "05d479f2b24b43debf234fc6ac9afa7c",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, description='OxidationStates', max=1181), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "training R2 = 0.989\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "728543ca398f40a5a4f72ff8e595003b",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, description='DensityFeatures', max=1181), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "59d108284fbc4116a5f420ba6b70a437",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, description='ElementProperty', max=1181), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "48b257ec6d594626a74e8ce94d42ee4f",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, description='OxidationStates', max=1181), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "training RMSE = 7.733\n"
     ]
    }
   ],
   "source": [
    "# make the pipeline\n",
    "pipeline = Pipeline([\n",
    "    ('featurize', featurizer),\n",
    "    ('regress', RandomForestRegressor(n_estimators=50, random_state=1)),\n",
    "])\n",
    "\n",
    "pipeline.fit(df, y)\n",
    "\n",
    "# get fit statistics\n",
    "print('training R2 = ' + str(round(pipeline.score(df, y), 3)))\n",
    "print('training RMSE = %.3f' % np.sqrt(mean_squared_error(y_true=y, y_pred=pipeline.predict(df))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cross Validation\n",
    "To run cross validation, the featurizer transformation can't be in a pipeline with the regressor, as the initial form of the data cannot be used with KFold. This is because the transformer adds and removes columns, it's more than just a simple function of the data. Instead the final featurized data can be computed beforehand, here as ``X``."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "542baff203b04e97b911b16561b88811",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, description='DensityFeatures', max=1181), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "07bf029044f546c797f19ed3ab2a6a6c",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, description='ElementProperty', max=1181), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "14fb18cb1caa4cafb0e2d96668275379",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, description='OxidationStates', max=1181), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "X = featurizer.transform(df)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define a KFold for cross validation. Using RepeatedKFold can reduce variance in the cross val score without increasing the number of folds, this is similar to bootstrapping, as the data is randomly subsampled multiple times by the KFold and then averaged. Using repeated folds is a good way to reduce variance if there is sufficient compute time to do so. For very computationally expensive models, such as DNNs, it is common to use a single train/validation split (not counting the excluded test data)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "crossvalidation = RepeatedKFold(n_splits=5, n_repeats=3, random_state=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is the same linear regression as before, now with RepeatedKFold cross validation. This gives a better estimate of how well the model will generalize than looking at training error. Cross validation usually gives a pessimistic estimate, and in practice the best performing model will be retrained on the full set of train/val data before testing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Cross-validation results:\n",
      "Folds: 15, mean R2: 0.899\n",
      "Folds: 15, mean RMSE: 22.894\n"
     ]
    }
   ],
   "source": [
    "lr = LinearRegression()\n",
    "\n",
    "scores = cross_val_score(lr, X, y, scoring='neg_mean_squared_error', cv=crossvalidation, n_jobs=1)\n",
    "rmse_scores = [np.sqrt(abs(s)) for s in scores]\n",
    "r2_scores = cross_val_score(lr, X, y, scoring='r2', cv=crossvalidation, n_jobs=1)\n",
    "\n",
    "print('Cross-validation results:')\n",
    "print('Folds: %i, mean R2: %.3f' % (len(scores), np.mean(np.abs(r2_scores))))\n",
    "print('Folds: %i, mean RMSE: %.3f' % (len(scores), np.mean(np.abs(rmse_scores))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is the same with the random forest regressor."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Cross-validation results:\n",
      "Folds: 15, mean R2: 0.924\n",
      "Folds: 15, mean RMSE: 19.652\n"
     ]
    }
   ],
   "source": [
    "# compute cross validation scores for random forest model\n",
    "rf = RandomForestRegressor(n_estimators=50, random_state=1)\n",
    "\n",
    "r2_scores = cross_val_score(rf, X, y, scoring='r2', cv=crossvalidation, n_jobs=1)\n",
    "scores = cross_val_score(rf, X, y, scoring='neg_mean_squared_error', cv=crossvalidation, n_jobs=1)\n",
    "rmse_scores = [np.sqrt(abs(s)) for s in scores]\n",
    "\n",
    "print('Cross-validation results:')\n",
    "print('Folds: %i, mean R2: %.3f' % (len(scores), np.mean(np.abs(r2_scores))))\n",
    "print('Folds: %i, mean RMSE: %.3f' % (len(scores), np.mean(np.abs(rmse_scores))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model Selection with Grid Search\n",
    "\n",
    "A pipeline can be used with Grid Search or Random Search for model selection and hyper-parameter optimization. This can include normalization, scaling, whitening, PCA / dimensionality reduction, basis expansion, or any other preprocessing or data transformation steps.\n",
    "\n",
    "Setting up a pipeline is design pattern that gives a straight forward abd repeatable method of processing the data and training a model. This can make it easy to try many different models and perform model selection with a hyper-parameter optimization scheme like grid search. \n",
    "\n",
    "Before doing model selection, the data should be split into a training set and a holdout test set. This tries to measure the generality of the model, predicting how it may perform on real data. Without a test set there is no way to measure if the model has likely overfit the training data.\n",
    "\n",
    "Note: The best model is chosen by cross validation score, and only the final model (after being retrained on all train/val data) is evaluated on the test set. Evaluating multiple models on the test set and choosing the best of them is an almost sure way of leading to overfitting or overestimating the true performance/generality of the model, exactly what we are trying to avoid by creating a hold out test set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9137335686081911\n",
      "{'n_estimators': 100}\n",
      "0.952679367459225\n"
     ]
    }
   ],
   "source": [
    "rf = RandomForestRegressor(n_estimators=50, random_state=1)\n",
    "param_grid = [\n",
    "  {'n_estimators': [10,15,20,25,30,50,100]},\n",
    "]\n",
    "gs = GridSearchCV(rf, param_grid, n_jobs=4, cv=5)\n",
    "gs.fit(X_train, y_train)\n",
    "print(gs.best_score_)\n",
    "print(gs.best_params_)\n",
    "print(gs.score(X_test, y_test))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Random Search\n",
    "Random search is another possible option for hyper-parameter selection, and usually outperforms grid search both in theory and in practice (see Random Search for Hyper-Parameter Optimization by Bergstra & Bengio). This is true especially in higher dimentional hyper-parameter spaces. This shows an example of a scaling step in the pipeline, which can improve performance for some types of models. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "best crossval score 0.901\n",
      "best params {'regress__n_estimators': 110}\n",
      "training R2 = 0.988\n",
      "training RMSE = 7.908\n",
      "test R2 = 0.953\n",
      "test RMSE = 16.281\n"
     ]
    }
   ],
   "source": [
    "pipe = Pipeline([\n",
    "    ('scale', StandardScaler()),\n",
    "    ('regress', RandomForestRegressor(random_state=1)), \n",
    "])\n",
    "\n",
    "param_dist = {'regress__n_estimators': sp_randint(10,150)}\n",
    "\n",
    "gs = RandomizedSearchCV(pipe, param_dist, cv=crossvalidation, n_jobs=-1)\n",
    "gs.fit(X_train, y_train)\n",
    "\n",
    "print('best crossval score ' + str(round(gs.best_score_, 3)))\n",
    "print('best params ' + str(gs.best_params_))\n",
    "\n",
    "# get fit statistics\n",
    "print('training R2 = ' + str(round(gs.score(X_train, y_train), 3)))\n",
    "print('training RMSE = %.3f' % np.sqrt(mean_squared_error(y_true=y_train, y_pred=gs.predict(X_train))))\n",
    "print('test R2 = ' + str(round(gs.score(X_test, y_test), 3)))\n",
    "print('test RMSE = %.3f' % np.sqrt(mean_squared_error(y_true=y_test, y_pred=gs.predict(X_test))))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "matminer python",
   "language": "python",
   "name": "matminer"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
